%
% Vasyl Mykhalchuk
% 14-02-12 
%
%
function [ strains ] = strainPortCPP( v1, v2, v3, v1_, v2_, v3_ )

% Middle of the v2-v3 edge 
center = (v2 + v3) / 2.0;
center_ = (v2_ +v3_) / 2.0;

% Gages vectors BEFORE the deformation
e_A = v2 - v1;
e_B = center - v1;
e_C = v3 - v1;

% Gages vectors AFTER the deformation
e_A_ = v2_ - v1_;
e_B_ = center_ - v1_;
e_C_ = v3_ - v1_;

% Compute gages e_p e_q e_r
e_p = ( norm(e_A_) - norm(e_A) ) / norm(e_A);
e_q = ( norm(e_B_) - norm(e_B) ) / norm(e_B);
e_r = ( norm(e_C_) - norm(e_C) ) / norm(e_C);

% Angles among gages: ang_q (btw. e_p and e_q) and ang_r (btw. e_p and e_r)
ne_A = e_A / norm(e_A);
ne_B = e_B / norm(e_B);
ne_C = e_C / norm(e_C);
cos_q = dot(ne_A, ne_B);
cos_r = dot(ne_B, ne_C);
ang_q = acos(cos_q); % Angle between e_A and e_B
ang_r = acos(cos_r); % Angle between e_B and e_C
ang_r = ang_q + ang_r;

%{
Normal strains (e_x, e_y) and shear strain (e_xy) are computed by solving a simple linear system.
%}
e_x = e_p;
% Matrices
left = zeros(2, 1); right = zeros(2, 2); unknowns = zeros(2, 1);
left(1, 1) = e_q - 0.5 * (1 + cos(2 * ang_q) ) * e_p;
left(2, 1) = e_r - 0.5 * (1 + cos(2 * ang_r) ) * e_p;
right(1, 1) = 0.5 * (1 - cos(2 * ang_q) );
right(1, 2) = 0.5 * sin(2 * ang_q);
right(2, 1) = 0.5 * (1 - cos(2 * ang_r) );
right(2, 2) = 0.5 * sin(2 * ang_r);
unknowns = right \ left;
e_y = unknowns(1, 1);
e_xy = unknowns(2, 1);

%{
Compute principal strains (e1 and e2) and the orientation from e_x, e_y, and e_xy above.
First, compute e_1, e_2, and phi (all scalars) using Mohr's circle
%}
R = 0.5 * sqrt( (2*e_xy) ^ 2 + (e_x - e_y) ^ 2 );
% !! assert(R>=0)
e_1 = ( e_x  + e_y) * 0.5 + R;
e_2 = (e_x + e_y) * 0.5 - R;

phi = atan( -(e_x - e_1) / e_xy );

n = cross(e_A, e_B);
n = n / norm(n);
e1 = RotVecArAxe(e_A, n, phi);
e1 = e1 / norm(e1);
e2 = cross(n, e1);

strains = [e_1, e_2];

end % function strainPortCPP

